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Abstract 



The problem of hydraulic fracture for the PKN model is considered within the framework of 
£Jjr^ approach presented recently by Linkov (2011). The modified formulation is further enhanced by 

employing an improved regularized boundary condition near the crack tip. This increases solution 
accuracy especially for singular leak-off regimes. A new dependent variable having clear physical 
sense is introduced. A comprehensive analysis of numerical algorithms based on various dependent 
variables is provided. 
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1 Introduction 

In its broadest definition, hydraulic fracturing refers to a problem of a fluid driven fracture propagating in 
a brittle medium. The process has been known for at least fifty years (Crittendon (1959), Desroches and 
Thiercelin (1993), Desroches et al. (1994), Detournay (2004), Geertsma and de Klerk (1969), Harrison et 
al. (1954), Hubbcrt and Willis (1957), Khristianovic and Zhcltov (1955)). The respective technology has 
been utilized in the petroleum industry to intensify the extraction of hydrocarbons for decades. Recently, 
due to economical reasons, it has been revived for exploiting non-conventional hydrocarbon deposits. 
The process has many other technological applications (e.g. disposal of waste drill cuttings underground 
(Moschovidis et al., 2000), geothermal reservoirs exploitation (Pine and Cundall, 1985) or in situ stress 
measurements (Desroches and Thiercelin, 1993). Hydrofracturing also appears in nature (e.g. geological 
processes, like magma-driven dykes (Rubin, 1995), (Lister, 1990) or a subglacial drainage of water (Tsai 
and Rice, 2010)). Due to the complexity of this multiphysical phenomenon, mathematical and numerical 
simulation of the process still represents a challenging task, in spite of the fact that immense progress 
has been made since the first algorithms were developed. 

The mathematical model of the problem should account for coupled mechanisms driving the process, 
which arc: i) solid mechanics equations, describing the deformation of the rock induced by the fluid 
pressure; ii) equations for the fluid flow within the fracture and the leak-off to the rock formation; iii) 
fracture mechanics equations, defining the conditions for fracture propagation. Further development of 
the model involves incorporation of mass transport for the proppant movement, fluid diffusion to account 
for the rock saturation by the leak-off flux, thermal effects affecting rhcological properties of the fluid 
and others. 
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The computational challenges of the hydrofracturing models result from several factors: (i) strong 
non-linearity introduced by the Poiseuille equation describing the fluid flow; (ii) in the general case, a 
non-local relationship between the fracture opening and the net fluid pressure; (iii) moving boundaries 
of the fluid front and the fracture contour (iv) possible lag between the crack tip and the fluid front. 

The first simplified mathematical descriptions of hydraulic fracture were summarized in the following 
three main classical models. The so-called PKN model was considered in Perkins and Kern (1961), 
where the authors adopted Sneddon's solution (Sneddon and Elliot, 1946) which was further enhanced 
by Nordgren (1972) to account for the fluid loss effect and fracture volume change. As a result, the crack 
length was determined as a part of the solution. The so-called KGD (plain strain) model was developed 
independently by Khristianovic and Zheltov (1955), and Geertsma and de Klerk (1969). Finally, the 
radial or penny-shaped model was introduced by Sneddon (1946) with constant fluid pressure and was 
extended for the general case by Spence and Sharp (1985). 

Different variations of the aforementioned models were used for treatment designs for decades, despite 
the fact that each of them is valid only under very specific assumptions (like elliptic cross-section of the 
fracture, and fracture half- length much greater than its constant height for the PKN model). Recently, 
the classical models have been largely replaced by the pseudo 3D models (Mack and Warpinski, 2000). 
A comprehensive review of the history and techniques of hydrofracturing simulation can be found in 
Adachi et al. (2007). 

Although the classical models have been superseded in most of the practical applications, they still 
play a crucial role when developing and analyzing new computational algorithms. The models enable 
one to understand the nature of, and reasons for, computational difficulties, find the remedies for them 
and extend these ideas to the general case. 

Thus, in the pioneering work by Nordgren (1972), main peculiarities of the model were analyzed and 
a numerical algorithm was proposed to deal with the problem. Asymptotic analysis of the solution near 
the crack tip for impermeable rock model was presented in Kemp (1989) and an approximate solution 
for the zero leak-off case, with an accuracy to the first four leading asymptotic terms, was constructed. 1 
Here, probably for the first time, the speed equation was efficiently implemented in the model. The 
fourth degree of the crack opening was considered as the proper variable and used in the numerical 
computations. Finally, the author showed that the numerical algorithm discussed in Nordgren (1972) 
corresponds to the Finite Volume (FV) scheme which was used in Kemp (1989). 

For the leak-off function defined by the Carter law the two leading terms of asymptotics can be found 
in Kovalyshen and Detournay (2009), where the PKN problem was revisited for the second time to take 
into account the multi-scale arguments in the spirit of Garagash et al. (2011). The authors also used 
the FV scheme to tackle the transient regime. 

Recently, the classical models of Nordgren and Spence & Sharp have been revisited again in Lin- 
kov (2011a) - Linkov (2011d). The author discovered that in some formulation the hydraulic fracture 
problem may exhibit ill-posed properties while the governing equation is singularly perturbed at the 
crack tip. To eliminate the difficulties resulting from these facts, a number of measures were proposed: 
i) speed equation to trace the crack front instead of the usually applied total flux balance condition 2 ; ii) 
the so-called e-rcgularization technique which consists of imposing the computational domain boundary 
at a small distance behind the crack tip; iii) new boundary condition to be imposed in the regularized 
formulation; iv) new dependent variables: the particle velocity and the crack opening taken in a degree 
to exploit the asymptotic behaviour of the solution; v) the spatial coordinates moving with the crack 
front and evaluation of the temporal derivative under fixed values of these coordinates. 

The advantages of the modified formulation were clearly demonstrated on the basis of developed 
analytical benchmarks in Linkov (201 Id). An immense improvement of solution accuracy, computation 

X A complete analytical solution for the impermeable rock case was constructed in terms of fast converging series in 
Linkov (2011c) 

2 Probably, for the first time, this idea was recalled in indirect way in Spence and Sharp (1985) and utilized by Kemp 
(1989), but later was abandoned as the particle velocity at the crack tip is difficult to compute numerically 
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efficiency and stability was shown. In Mishuris et al. (2012) a further step in employing the modified 
formulation was done by analyzing the stiffness of the system of differential equations arising after 
spatial discretization. An efficient modification of the algorithm has been proposed to trace the fracture 
propagation. Finally, the investigation of solution sensitivity to some technical parameters has been 
considered. 

However, the aforementioned analysis is concerned with the situations when the leak-off vanishes 
near the crack tip and the fluid velocity does not change much along the entire fracture. 

The primary aim of this paper is to verify the recipes delivered in Linkov (2011a, b,c,d) and Mishu- 
ris et al. (2012) for an arbitrary leak off regime. In particular, the analysis includes: (i) Investigation 
of performance of numerical algorithms for hydrofracturing based on different formulations (different 
dependent and independent variables) ; (ii) Utilization of a new dependent variable defined as an integral 
of the crack opening. Such a variable has a clear physical and technological sense and is not related 
to any specific type of solution asymptotic behaviour near the crack tip; (iii) Modification of the way 
to impose boundary conditions in the framework of the e-regularization technique; (iv) Identification 
of optimal ranges of the various technique parameters, accuracy of computations and stiffness of the 
resulting dynamical systems. 

The structure of this paper is as follows. In the next section we collect known results for the PKN 
model in various formulations. We restrict ourselves only to the information which is absolutely necessary 
to understand the paper. Subsection 2.3.2 contains an alternative formulation in terms of a new proper 
dependent variable - fracture volume. Also, some new results concerning the asymptotic behaviour of 
the solutions for different leak-off regimes are presented in Appendix A and Appendix B. Finally, since 
different authors use various notations, we adopt a unified style in our paper. 

In Section 3 we discuss in details the numerical procedures and present main results of the computa- 
tions comparing performances of different solvers under considerations. Finally, the main findings of this 
paper are summarized in the conclusion section. All these amendments are presented on the example of 
the PKN model but could be useful for other models of hydraulic fracture. 

2 Problem formulation and preliminary results 

2.1 Physical fundamentals and basic equations 

Consider a symmetrical crack of the length 21 situated in the plane x E (— i, I), where the length I = l(t) 
is one of the solution components changing as a result of the fluid flow inside the crack. The initial 
crack length is assumed to be nonzero: 1(0) = Z* > 0. There are reasonable motivations behind this 
assumption. Namely, for the initial (unstable) stage of the crack propagation the acceleration of the 
process is too large to neglect the inertial terms. For this reason, any classical model not accounting for 
this effect is not credible. On the other hand, in many cases the hydrofracturing process is associated 
with the so-called 'perforation technique'. The latter consists of the creation of a number of finger- 
shaped initial fractures by detonations of shaped charges spaced along the wellbore (Economides, 2000). 
In this way, hydrofracturing starts simultaneously from a number of non-zero length cracks. Moreover, 
in many rock formations (e.g. shale reservoirs) cracks already exist but are closed by the confining stress 
(Economides, 2000). Finally, the uncertainties involved in this complex multiphysics problem itself do 
not allow one to make any reliable modelling of the crack nucleation in the rock formation. 

By convention, we assume that the crack is fully filled by a Newtonian liquid injected at known rate 
qo(t) at the crack mouth x = 0. 

The Poiseulle equation for the Newtonian liquid flow in a narrow channel is written in the form: 



where 



where q = q(t, x) is the fluid flow rate and w = w(t, x) is the crack opening, while p = p(t, x) is the net 
fluid pressure, that is, the difference between the fluid pressure pj inside the fracture and the confining 
stress uo (p = pf — (To). The constant M involved in the equation is defined as M = I2(i, where [i stands 
for the dynamic viscosity (see for example Linkov (2011d)). 

The continuity equation, accounting for the crack expansion and the leak-off of the fluid, may be 
expressed as: 

^ + |^+«=0, t>0, 0<x<l(t), (2) 

where qi = qi{t,x) is the volume rate of fluid loss to formation in the direction perpendicular to the 
crack surfaces per unit length of fracture. 

Numerical algorithms for the PKN model, with and without leak-off, have been considered in Nord- 
gren (1972), Kemp (1989), Kovalyshen (2010) and others. Improvements basing on the speed equation 
and the e-regularisation technique have been introduced in Linkov (2011c, d); Mishuris et al. (2012), for 
the case when leak-off vanishes near the crack tip. Below, we discuss the effectiveness of this approach on 
the three most popular leak-off models: Carter law (Carter, 1957), modified law incorporating pressure 
difference (Clifton and Wang, 1988) and bounded leak-off near the crack tip. An extensive discussion 
on possible behaviour of the leak-off function can be found in Kovalyshen (2010). 

In our numerical simulations, we utilise one of the following leak-off variants: 

qi{t,x)=q ( f\t,x)+q*{t,x), j = 1,2,3, (3) 

^ (1) = ^=4v = ,f = c^{t) P+ c^t), o< *<«*)• (4) 

V* — T \ x ) yt — T{x) 

Here C\ = Cl is usually assumed to be a known constant defined experimentally (Carter, 1957). Recently 
it was estimated analytically for a poro-elastic material in Kovalyshen and Detournay (2009). The 
function t(x) contains information on the history of the process. It defines the time at which the 
fracture tip reaches the point x and can be computed as the inverse of the crack length: 

t(x) = r x {x), x>h. (5) 

For x < U we conventionally set t(x) = 0. Other constants in (4), Cj(t) = Cj(t,w,p), (j = 2,3) may 
depend on the solution itself but arc bounded functions in time. Finally, we assume that the terms q* , 

(j = 1,2,3) in (4) are negligible in comparison with q\^ near the crack tip. Note that application of 
the Carter leak-off law (Carter, 1957) which is a simplified model of established fluid diffusion through 
the fracture walls, is rather questionable at the initial stage of the crack propagation (Nordgren, 1972; 
Lenoach, 1995; Mathias et al, 2009). 

In this paper, we are aiming to build a numerical scheme in a general form regardless of the particular 
choice of leak-off function. Thus, the collection of possible leak-off representations given in (4) covers 
the whole spectrum of possible bahaviours used in the hydrofracturing simulations (Kovalyshen, 2010). 

The system of equations (1) - (2) should be supplemented by the elasticity equation. We consider 
the simplest relationship used in the PKN formulation 

p = kw, (6) 

with a known proportionality coefficient k — ^ -i found from the solution of a plane strain elasticity 
problem for an elliptical crack of height h (Nordgren, 1972). Constants E and v are the Young modulus 
and the Poisson ratio, respectively. In fact, (6) implies that the so-called zero toughness condition is 
satisfied in this model (Kic = 0). 
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On substitution of the Poisculle equation (1) and elasticity relationship (6) into the continuity equa- 
tion (2), one obtains a well known lubrication (Reynolds) equation defined in the trapezoidal domain 
(t > 0, < x < l(t)): 

dw k d ( % dw\ 

w 3 —)+qi=0. (7) 



dt Mdx V dx. 

Since the system has its natural symmetry with respect to variable x and the equations are local, 
it is convenient to consider only half (symmetrical part) of the interval [0, l(t)] instead of the full crack 
length [-l(t),l(t)]. 

Following the discussion on the initial crack length above, the initial conditions for the problem are: 

1(0) = h, w(0,x) = w*(x), x€(0,U). (8) 

The boundary conditions include: i) known fluid injection rate at the crack mouth, h) zero crack 
opening and zero fluid flux rate at the crack tip: 

q(t,0) = q Q (t), w(t,l(t))=0, q(t,l(t))=0. (9) 

Note that the problem formulated in this way looks overdetermined as the governing equation (7) is of 
the second order with respect to spatial variable. This issue shall be discussed later. 

Finally, by consecutive integration of equation (7) over time and then space, one can also derive the 
standard formula for the global fluid balance in the form: 

l(t) ,t Mt) ,t 

[w(t,x) - w*(x)}dx - / q (t)dt+ / / qi(t, x)dtdx = 0, (10) 



where it is accepted that w*(x) = when x > i* and l'(t) > 0. 

As has been shown in Linkov (2011d), the crucial role in the analysis of the problem plays the particle 
velocity defined in the following manner: 

V(t,x) = -, i>0, 0<x<l(t), (11) 
w 

which indicates the average velocity of fluid flow through the cross-sections of the fracture. 

Under the assumption that the crack is fully filled by the fluid and sucking, ejection or discharge 
through the front can be neglected, the fluid velocity defines the crack propagation speed and the 
following speed equation is valid (Kemp, 1989; Linkov, 2011a, b) 3 : 

l'(t) = V(t,l(t)), t>0. (12) 

Moreover, for physical reasons, one can assume that the fluid velocity at the crack tip is finite 

< V(t,x) < co, t>0, x<l(t) (13) 

Note that, allowing the crack propagation speed to be infinite, one has to simultaneously include the in- 
ertia term in the equations. Thus, the estimate (13) is a direct consequence of neglecting the acceleration 
terms. 



3 In fact, the speed equation in this form is valid only under the assumption of zero spurt loss at the crack tip (Nordgren, 
1972; Clifton and Wang, 1988; Adachi et al., 2007) 
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2.2 Asymptotic behaviour of the solution and its consequences 

As was mentioned in Spcncc and Sharp (1985), the fact that both w and q are present in (11), creates 
serious difficulties when trying to use the fluid velocity as a variable. However, as was shown in Linkov 
(2011a,b,c,d), proper usage of fluid velocity may be extremely beneficial. First, it allows to substitute 
two boundary conditions at the crack tip (9)2,3 with a single one additionally incorporating information 
from the speed equation (12), (13). 

Indeed, the boundary conditions (9)2,3 in view of (1) and (6) lead to the estimate 

w(t,x)=o((l(t)-x)i), x-*l(t), (14) 

which does not necessarily guarantee (13). However, further analysis of the problem, for different leak-off 
functions (see Kemp (1989); Kovalyshcn and Detournay (2009) and Appendix B of this paper), shows 
that the particle velocity is bounded near the crack tip and the crack opening exhibits the following 
asymptotic behaviour: 

w(t, x) = wo(t)(l(t) - x) 3 + Wl (t)(l(t) - x)° + o((/(i)-x) Q ), as x^l(t), (15) 

with some a > 1/3. For the classical PKN model for an impermeable solid (or when leak-off vanishes 
near the crack tip faster than the crack opening) the exponent a = 4/3 was found in Kemp (1989) and 
Mishuris et al. (2012). For the case of the singular Carter's type leak-off, the exponent a = 1/2 was 
determined in Kovalyshen and Detournay (2009). 

Note that the asymptotics (15) shows that fluid velocity is indeed bounded near the crack tip. 
Moreover, 

V(t, x) = V (t) + Ki(t)(i(t) -xf + o((l(t) - xf), (16) 
as x — > l(t), where (3 = a — 1/3 and 

V o = ^o(t), V 1 = ^{a+^\wl{t)w 1 {t). (17) 

As follows from Appendix B, V(t,x) may not be so smooth near the crack tip as one could expect and 
the exponent (3 in (16) plays an important role for this. Indeed, if (3 > 1 then V(t, •) £ C^O, L(t)] and the 
particle velocity function is smooth enough near the crack tip. However, this happens only in the special 
case of a = 4/3 when V x (t,x) is bounded near the crack tip. In case of singular leak-off (0 < j3 < 1), 
the particle velocity near the crack tip is only of the Holder type V(t,-) G C* 1 [0, L(t)) f] H^[0, L(t)]. In 
Appendix B we present an exact form of the asymptotic expansion (15), which yields the aforementioned 
smoothness deterioration of V near the crack tip for the singular leak-off models. 

Note that estimate (15) (or equivalently (16)) enforces the condition (13). Thus, in view of (11), the 
pair of conditions (9)2 and (9)3 is equivalent to (9)2 and (15). This discussion clearly illustrates why 
accounting for asymptotic behaviour of the solution in form (15) is of crucial importance for effective 
numerical realisation of any algorithm utilised in hydrofracturing. On the other hand, the fact that the 
particle velocity function is smooth enough near the crack tip has been one of the important arguments 
to use the speed equation and proper variable approach as the basis for improvement of the existing 
numerical algorithms (Linkov, 2011a). It should be emphasized that behaviour of V(t,x) near the crack 
tip may have serious implications when using e-regularization technique. 

Therefore, the main aim of this paper is to show that, regardless of possible smoothness of the particle 
velocity near the crack tip, the approach proposed in Linkov (201 Id); Mishuris et al. (2012) is still 
efficient. 
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2.3 Normalised formulation 

Let us normalize the problem by introducing the following dimensionlcss variables: 

*ntfo(*) = ?o(*), i n qi(t,x) = Uqi(t,x), (18) 
where a; G (0, 1) and £(0) = 1. Using this notation, one defines the normalised particle velocity as: 

u> 2 <9?7j 



The Reynolds equation (2) in the normalised domain is rewritten in the following manner: 



dw 



(xV{t,l)-V{i,xj) 



dw _ 9V" 
dx dx 



Qi(t,x), 



at L(t) 

The leading terms of the asymptotic estimate of the leak-off function from (4) are now: 



(19) 



(20) 



Vl — x VI — £ 



g, (3) (t,i) - C- 31 (*>(£,£) + C 32 (t). 



Here, the function 



D(t) = 



L(t)> 



(21) 



(22) 



is introduced in the Appendix A, where the remainder between the normalised total flux and the leading 
term (21) has been effectively estimated. Thus the normalised term q*(t,x) vanishes near the crack tip 
faster than the solution itself. 

Finally, normalised initial conditions (8) and boundary conditions (9) are: 



and 



L(0) = 1, w(0,x) = w*(x), a: €(0,1), 



-j^ r) w 3 ^(t,0) = q (t), w(t,l) = 0. 



The global fluid balance (10) can be written in the form: 



L(t) / w(t,x)dx 



{x,0)dx- / q (t)dt+ / L(t) / qi(t,x)dxdt = 



(23) 
(24) 

(25) 



For convenience, from this point on we will omit the "~" symbol for all dependent and independent 
variables and will only consider the respective dimensionless values. 

Note that the particular representation (20) of the Reynolds equation highlights an essential feature 
of the problem - it is singularly perturbed near the crack tip. Indeed, both coefficients in front of the 
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spatial derivatives on the right-hand side of the equation (20) tend to zero at x = 1. Thus, the asymptotic 
behaviour of the solution near the crack tip (x — > 1) 



w 



w (t) (1 - x)* + Wl (t) (1 - x) a + o ((1 - x) a ) , 

V = V (t) + Vt(t) (1 - x) a ~^ + o ((1 - a:)«-*) , 
represents nothing but the boundary layer. Moreover, normalizing (17)i one obtains: 

V (t) 



3L(t) 



w 3 (t). 



(26) 
(27) 

(28) 



Note that terms wo, wi and Vo, V\ in (26) and (27) are different from those in (15) and (16). In fact, 
the former should be written with "~" symbol. 

On substitution of (19) into (20), one eliminates the particle velocity function from Reynolds equa- 
tion: 



dw 
~dt 



1 



L 2 (t) 



1 3 dw , o~2 

-w x— + 3w 



dw 
dx 



,d 2 w 
dx 2 



ox- 



(29) 



Here wq is the multiplier of the first term of the asymptotic expansion (15). This form of lubrication 
equation exhibits the same degenerative properties as (20). Also the coefficients appearing in front of 
the leading terms tend to zero near the crack tip. 

The speed equation (12) defining the crack propagation speed is given in the normalised variables as: 



L'{t) = V (t), t > 0. 
Taking into account (28), the latter can be rewritten in the following form 



t > o. 



(30) 



(31) 



This equation serves us to determine the unknown value of the crack length L(t). As it has been shown 
in Mishuris et al. (2012), such an approach has clear advantages over the standard one based on the 
global fluid balance equation (25). 

As a result of the foregoing transformations, one can formulate a system of PDEs describing the 
hydrofracturing process. The system is composed of two operators: 



dt 



w = A w (w,L 2 ), -L 2 {t)=B w (w) 
at 



(32) 



where A w is defined by the right-hand side of equation (29) with the boundary conditions (24)1,2, while 
the second operator B w is given by (31). The system is equipped with the initial conditions: 



L(0) = 1, w(0,x)=w*(x), a; €(0,1). 



(33) 



2.3.1 Reformulation of the problem in proper dependent variable. First approach 

In Linkov (2011d) and later in Mishuris et al. (2012) it has been shown that the dependent variable 

U(t,x) = w 3 (t,x) (34) 

is more favorable for the solution of the system (32), (33) than the crack opening itself. This idea is 
based on the fact that, according to the asymptotics of the solution near the crack tip, the dependent 
variable U is much smoother than w. In the case of an impermeable solid, the solution U is analytic 
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in the closed interval [0, 1] (see Linkov (2011d)). However, the type of leak-off function is of significant 
importance here. Thus, adopting asymptotic representation (26), one can see that for x — > 1 



U = U (t)(l - x) + U x {t){l - x)i +a + o((l - x)i +2a ), 



(35) 



where the coefficients Uo(t) and U\(t) are directly related to those appearing in the crack opening 
formulation: 

U (t) = w 3 (t), Ut(t) = w 2 (t) Wl (t). (36) 

Depending on the type of leak-off described in (4), the exponent in the second asymptotic term | + a 
may take value 3/2, 11/6 or 2, respectively. Thus in the first two cases, the transformation (34) no longer 
results in polynomial representations of asymptotic expansion for U. For this reason, the advantage of 
the approach using variable U in more general cases, when the leak-off is singular near the crack tip, 
should still be confirmed. This is one of the aims of this paper. On the other hand, at least two factors 
work in favor of this formulation. First the first derivative of U is not singular and, second, the particle 
velocity is given by a linear relationship 



V(t,x) 



i du 

3L(t) dx ' 



(37) 



The governing equation (20) in terms of the new variable can be written in the normalized domain 
x e (0, 1) as: 

dU _ _L_ 
~dt ~ 



L(t) 



dU dV 
(xV(t,l)-V(t,x))^--3U^- 



(38) 



Similarly to (29) the particle velocity function may be eliminated from the lubrication equation: 



8U_ 
~dt 



3L 2 (t) 



Xllr, 



dU 
dx 



dU 
dx 



3U 



,d 2 U 
dx 2 



-3U*q t . 



(39) 



Note that equations (38)-(39) are of a very similar structure to those evaluated for the crack opening w. 
They exhibit the same degenerative nature near the crack tip. 
Finally, boundary conditions (24) transform to: 



3L(t) dx 

while the speed equation (31) takes the following form 

2 



U(t,0) = q (t), t/(M)=0, 



—L 2 

dt 



;Uo(t), 



t > 0. 



The system of PDEs equivalent to (32) is now defined as: 



dt 



U = Au(U,L 2 ), 



dt 



L 2 (t)=Bu(U). 



(40) 



(41) 



(42) 



The operator An is described by (39) with boundary conditions (40), while the second operator Bu is 
given by (41). Finally, the initial conditions are similar to those in the previous formulation (23): 



L(0) = 1, U(0,x)=wl(x), x 6(0,1). 



(43) 
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2.3.2 Reformulation of the problem in proper dependent variable. Second approach 

The aforementioned formulation of the problem in terms of the dependent variable U has one considerable 
drawback. It is well known that for different elasticity models and different hydrofracturing regimes one 
has various asymptotic behaviours of the solution near the crack tip (Adachi and Dctournay, 2002). 
For example, for exact equations of elasticity theory and the zero toughness condition (Kic = 0), the 
exponent of the leading term of if varies from 2/3, for the Newtonian fluid, to 1, for the ideally plastic 
fluid. Thus, the same reformulation to the type of the proper variable might be inconvenient, or even 
impossible. 

For this reason, we introduce another dependent variable. Although it does not transform the 
asymptotic behaviour of the solution in such a smooth manner as it has been done previously when 
adopting U, this variable has its own advantages. Namely, let us consider a new dependent variable 
defined as follows: 4 

Sl(t,x)= [ w(t,S)dt. (44) 



This variable is not directly related to any particular asymptotic representation ofw(x,t), however it 
assumes that w — > for x — ^ 1. As a result the form of governing equations for fl remain the same 
regardless of w(x,t) asymptotics, i.e. this formulation has a general (universal) character. Note, that 
in case of U the optimal way of transformation for the lubrication equation essentially depends on the 
exact form of asymptotic expansion (the leading term) for w. 

Another advantage of Q comes from the fact that it has clear physical and technological interpretation. 
Namely, it reflects the crack volume measured from the crack tip. 

Asymptotics of the function O near the crack tip takes the following form 5 : 

n(t,x) = n (t)(i-x)^ +n 1 (t)(i-x) a+1 + o((i-x) a+1 ), x->i, (45) 

where the coefficients f2o(*) an d f2i(t) are related to those in (26): 

fi (t) = ^w (t), Qi(t) = ——wi(t). (46) 
4 a + 1 

Thus, similarly to U , the new variable is smoother than the crack opening, w, and the first singular 
derivative of Q is that of the second order. 

By spatial integration of (20) from x to 1 and substitution of (44) one obtains: 



dfl 1 



dt L(t) 



(v(t,x)-xv(t,i)) — + v(t,i)n 



Qi, (47) 



where the monotonicity of L'(t) > has been taken into account and 

Ql(t,x) = f ®(t,£K- 

J X 

Here, the particle velocity (19) is computed in the manner: 



VM-^l(£)- («> 



3L(t) dx \ dx 

By eliminating V(x, t) from the equation (47) we derive a new formula for the lubrication equation 



dt L 2 (t) 



on\ 3 & 2 n 64 3 / an 

dx~) ^ + si "r^fe 



Qi- (49) 



4 This variable is compatible with the idea of the FV scheme. 

5 Interestingly, x) behaves near the crack tip (x — * 1) similarly to w 4 (t, x) used by Kemp (1989). 
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The boundary conditions (24) are expressed in the following way: 



^VS«,0) = ,„, g(M) = 0. ,50, 



L(t) \ dx J dx 2 ' ' dx 

Interestingly, the first boundary condition, directly substituted into the lubrication equation (47) can be 
equivalently rewritten in the form 

§ft°>--5^^raw+t$-<j*o>, (51) 

This condition, in turn, represents nothing but the local (in time) flux balance condition. To verify 
this, it is enough to apply the time derivative to the equation (25). Furthermore it appears much easier 
for implementation into a numerical procedure than (50) i itself, but may lead to some increase of the 
problem stiffness, as we will show later. 

It is easy to check, by using the governing equation (47) and limiting values of all its terms for x — > 1 
, that a weaker boundary condition 

Q(M) = (52) 

is equivalent to the original one (50)2- Finally, the speed equation (31) in the formulation assumes 
the following form: 

- (53) 

In this way we obtain another system of PDEs that is composed of two operator relations: 

jn^Ani^L 2 ), j t L 2 (t)=Bn(n), (54) 

where, as previously, Aq is defined by (29) with boundary conditions (50)1,2 or (51) and (52). The 
second operator, Bq, is given by equation (53). Here the initial conditions are obtained from (23): 

L(0) = 1, n(0,x) = Sl*(x) = I tu*(f)eZ£. (55) 

J X 

2.4 -s-regularization and the respective boundary conditions 

In our analysis we are going to use the so-called e-rcgularization technique. It was originally introduced 
in Linkov (2011a) for the system of spatial coordinates moving with the fracture front. In Mishuris et 
al. (2012), the authors efficiently adopted the approach for the normalised coordinate system. 

The reason to separate the domain from the end point x = 1 by a small distance of e and to introduce 
e-regularisation has been thoroughly described in Linkov (2011a). It consists of replacing the Dirichlet 
boundary condition (40)2 with an approximate one: 

U[t, 1 -e) = 3eL(t)V(t, 1), (56) 

emerging from deep physical arguments. The value of the crack propagation speed V(t, 1) (and simul- 
taneously the particle velocity at a fracture tip) was suggested to be computed from the speed equation 
(41) in its approximated form: 

^^-m^' 1 -^ (57) 

The pair of conditions (56) - (57) has shown an excellent performance in terms of solution accuracy and, 
as has been proven in Mishuris et al. (2012), reduced the stiffness of dynamic system in case of leak-off 
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function vanishing near the crack tip. One can check that for such a leak-off model numerical error 
introduced by using the approximate conditions, instead of the exact ones, is of the order 0(e 2 ). In view 
of all improvements following from utilisation of the regularized conditions (Linkov (2011d); Mishuris et 
al. (2012)) such a strategy is fully justified and in fact inevitable. 

The conditions can be written in an equivalent form. Indeed, one can merge (56) and (57) into a 
single condition of the third type: 

dU 

U{t,l-e )+£— (t,l-e)=0. (58) 
Ox 

Interestingly, the latter condition is nothing but the consequence of a direct utilization of the information 
about the leading term of asymptotics of the solution near the crack tip (compare with (35)). 

Analogously, one can define the respective pairs of boundary conditions in the regularized formula- 
tions. Considering the dependent variable uu one should take (31) together with the condition 

w(t,l-E)+3e-^-(t,l-e) = 0, (59) 
Ox 

while analysing the system based on the dependent variable 57, the speed equation (53) should be 
accompanied by 

4Sl(t, l-e) + 3e— (t, 1-s) = 0. (60) 
ox 

To conclude this subsection, one can make a prediction that in the case of singular leak-off function, 
even when using the e-regularization technique, the accuracy of solution should be worse than that 
presented in Linkov (201 la), Mishuris et al. (2012). However, it is always possible to use information on 
accurate asymptotic behaviour of the solution (employing higher order terms) and in this way improve 
the accuracy of computations. 



3 Numerical solution of the dynamic systems 

In this section, three alternative systems of PDEs (32) , (42) or (54) describing the problem of hydrofracru- 
ing are transformed into the corresponding non-linear dynamic systems of the first order. Then, on the 
basis of respective analytical benchmarks, we analyze their stiffness properties, the accuracy and effi- 
ciency of computations. The benchmark solutions in question are described in Appendix C. 

3.1 Representation of the boundary conditions and the speed equation 

Consider a spatial domain of the problem reduced in accordance with the e-regularization technique to 
the interval x E [0, 1 — e], where e is a small parameter. Let the mesh, {xj}^ =11 be composed of N nodes 
with x\ = and xn = 1 — e. 

For each of the problem formulations, two boundary conditions should be accounted for: one specified 
at the crack inlet and a regularized boundary condition at x = 1 — e. In the following we present a brief 
description of how these conditions are introduced to the numerical scheme. 

From now on, for the dependent variables discussed above (w(t,x), U(t,x), Q(t, x)), we use common 
notation f(t,x) together with a convention = f(t,Xk)- 

To discretize the first boundary condition (depending on the formulation: (24)i, (40)i or (50)i) 
we exploit the smooth character of the solution near the point x = 0. Thus, accepting a polynomial 
approximation of f(x, t) on the interval x € [x±, X3], the respective nonlinear relation between /1, f 2 and 
/3 may be derived: 

A 1 (f 1 ,t)f 1 + A 2 (f 1 ,t)f 2 + A 3 (f 1 ,t)h = q a . (61) 
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Naturally, one can use a two-point FD representation of the condition at x = 0. However, for consistency 
and a better numerical performance, formula (61) is consequently used in this paper. The only exception 
is made in the case where the boundary condition in formulation (51) is utilized. The latter shall be 
discussed later on. 

As mentioned in 2.4, the regularized boundary condition in the e-regularization technique proposed in 
Linkov (2011a) is equivalent to a mixed boundary condition based on the leading term of the asymptotic 
expansion (see (58), (59), (60)). Below we propose a modification of this approach which consists in 
employing two terms of the asymptotics. We will show that such strategy prevents the deterioration 
of accuracy when the solution is not so smooth as in the cases originally considered in Linkov (2011a), 
Mishuris et al. (2012). 

According to (26), (35) and (45), the following asymptotics approximation 6 is acceptable in the 
proximity of the crack tip (x G [xn-2, 1]): 

f(t,x) = e[ f \t)(l - xr + e{ f \t)(l - x) a \ (62) 

The values of ol\ and a 2 are known in advance and depend, as it has been discussed above, on the chosen 
variable and the behavior of the leak-off function. Assuming that the last three points of the discrete 
solution (xat_2, /at_2), (xjv-i, /iv-i) and (xn,In) lie on the solution graph (x, f(t, x)), one can derive 
a formula combining all these values in one equation: 

In + b^ljN-i + = 0, (63) 

where b- = bj (xn-2, xn-i, xn)- Relation (63) is consequently used to represent the regularised 
boundary condition at x = 1 — e. 

Remark 1. In the authors opinion, the presented approach is a direct generalization of that proposed 
in Linkov (2011a). Indeed, if one takes e 2 = in the representation (62) then the pair of the equations 
(58) and (57) follows immediately. If a 2 — a\ = 1, which means that the leak-off function qi is 
bounded near the crack tip, the second asymptotic term provides a small correction. However, in the 
case of the Carter law, when a 2 ~ a i = 1/6, it brings an important contribution and improves the 
accuracy of the computations, as will be shown later. 

Finally, coefficient e\ from (62) is substituted into the pertinent form of the speed equation (31), 
(41) or (53) to give the ordinary differential equations for the crack length: 

Note that the right-hand sides of the equations define the boundary operators B w , Bjj and 2?n from 
(32) 2 , (42) 2 or (54) 2 , respectively. 

Remark 2. As it follows from this analysis, the e-regularization is, in a sense, equivalent to the 
introduction of a special tip clement in the discrete solution. Thus, one can see a complementarity with 
the approach utilised in Kovalyshen (2010). However, and this is crucial for the analysis, only the speed 
equation together with e-regularization allows to take into account both the local and global phenomena, 
and to do this in the most efficient way from the numerical point of view. 

Remark 3. In the case of the dependent variable U, apart from the representations (62) of the 
boundary condition near the crack tip in the linear form 

U N = b<f J) U N - 1 + b¥ ) U N - 2 , (65) 

one can use a nonlinear one, adopting the relationship between this dependent variable and the crack 
opening w: 

U N = (b[ w) + b[ w) i/u^~ 2 ) 3 . (66) 

6 One can argue that (62) is used as a test function for the crack tip element. 
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The fact that the latter condition is nonlinear is not essential for the implementation, as the entire 
problem (boundary condition at the point x = and the governing equation itself) is nonlinear. On 
the other hand, the two terms representation (65) of the function U is less informative than the same 
representation for the functions w (or Q) and thus, using the modified condition (66), one can expect a 
better solver performance. 

3.2 Spatial discretization of the Reynolds equation. Corresponding dynamic 
systems 

Let us consider the Reynolds equations written in different dependent variables ((32), (42) or (54), re- 
spectively) . By representing the spatial derivatives in the right-hand sides of the corresponding equation 
by central three point finite difference schemes, we obtain a nonlinear system of TV — 2 ordinary dif- 
ferential equations for the values fi(t) at each internal point of the spatial domain (x2, a;jv-i). The 
respective boundary conditions are embedded into the system through equations (61) and (63). 

Supplementing the system with the pertinent form of the speed equation (64), we obtain a non-linear 
dynamic system of first order describing the process of hydrofracturing which can be written in the form: 

F'=A<^F + G<fl, (67) 

where F = F(t) is a vector of unknown solution [f(t, x\), f(t, X2), ■ ■ ■ , f(t, xn), L 2 (t)] of dimension TV — 1. 
Note that matrix A^) and vector G^- 1 depend, generally speaking, on the solution. Matrix A is the 
so-called mass matrix of the system, in the case of which a tri-diagonal form prevails (however the 
boundary conditions and the last equation for L 2 (t) disturb the tri-diagonal structure). 

Remark 4 In the case when the boundary condition in formulation (51) is in use, the dimension of 
the dynamic system is TV. Indeed, this condition has the form of an ODE, and thus can be substituted 
directly into the system as an additional equation. 

In our numerical computations two different types of spatial meshes are used. The first one is a 
regular mesh (with uniformly spaced nodes): 

a$ = (l-e)^, m=l,...,TV. (68) 

The second mesh gives an increased node density when approaching the crack tip (8 > 1): 

= l-(l-(l- £ *)^)', m = l,...,TV. (69) 

In the foregoing index m refers to the number of the nodal point. Mesh (69) allows one to choose 
appropriate parameter 8 to suppress the stiffness of dynamic system or to increase the solution accuracy. 

The stiffness of a dynamic system may be described by the condition number or the condition ratio 
(Aiken, 1985) of a mass matrix A^. In general, the values given by various measures are different (see 
some consequences in Mishuris et al. (2012)). In this paper we use the condition ratio as the measure of 
the system stiffness. Some rough estimation of this parameter may help to choose an optimal value of 8 
from the stiffness point of view. In the case under consideration, one can analyse the condition ratio of 
a simplified variant of the system (39), where only the leading term (with the second order derivative) is 
preserved and the nonlinear multiplier is substituted by the first term of the asymptotic expansion for 
U. It turns out that 6 = 2 gives the lowest possible stiffness. 7 Naturally, in the general case, the optimal 
value of 5 can be different. We have checked however that, for three alternative problem formulations, 
there are three different optimal values of the parameter, but each of them is very close to 2. Thus, in 
the following section all results concerning nonuniform mesh are presented for 6 = 2. 

interestingly, the same exponent <5 = 2 was used due to completely different reason in Spcnce and Sharp (1985) when 
the cost functional was minimized on the entire domain (0,1) (thus, e = 0). 
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a) b) c) 




N N N 



Figure 1: Stiffness ratio of the matrix from (67) presented as a function of the number of the 

nodal points N in logarithmic scale for various dependent variables (w, U and Q): a) Benchmark cp 1 ) 
(Qi/lo — 0.9) for the uniform mesh, b) Benchmark q^ 1 ' (Qi/qo = 0.5) for the same mesh and and c) 
The benchmark (Qi/qa = 0.5) for nonuniform mesh 

s^rn with 5 — 2. For all cases we use the same 
regularised distance from the crack tip defined by e = 10~ 3 . The only difference between the systems 
for f2(i) and f2( 2 ) is the utilised boundary condition at the zero point: (51) and 61), respectively. 



3.3 Stiffness analysis 

In our analysis, we quantify the system stiffness by a condition ratio ka defined in (70). Since the 
problem is nonlinear, our investigation is to be done for the linearized form of matrix It is obvious 

that for all six variants of benchmark solutions under consideration (see Appendix C) one has different 
values of A^'. Computations are carried out for two types of meshes (the uniform and non-uniform 
one). For each of the benchmark solutions one obtains a constant (independent on time) value for the 
condition ratio apart from the fact that the matrix depends on time. Those values of the condition 
ratio ka are, generally speaking, different for various benchmarks and chosen meshes. 

Before comparing the results for various dependent variables, we have checked that for the dynamic 
system based on U the stiffness is almost identical for both forms of the regularized boundary condition 
at x = 1 — e ((58) and (63) respectively). Thus, for the rest of the dependent variables (w and fl) 
we restrict our stiffness investigation only to the formulation (63). Remarkably, the situation changes 
dramatically when one considers the accuracy of computations, which shall be discussed later on. 

Computing the condition ratio for the next variants of the matrix linearization, we have found out 
that the following estimate for large value of A" is valid for all cases under investigation: 

= ^4 ~ zu^N 2 , N->oo. (70) 

Here |A maa; |, and |A m i„| are the largest and smallest absolute values among the A*-^ matrix eigenvalues, 
while the constant vu^' is to be estimated numerically. Its values for all six benchmark cases are shown 
in Table 1. 

The representation (70) has been analytically derived in Mishuris et al. (2012) for the systems based 
on the dependent variable U . 

The following analysis includes investigation of stiffness sensitivity to: i) the solution (benchmark) 
type, ii) choice of the dependent variable, iii) choice of the independent variable (spatial mesh), iv) value 
of the rcgularization parameter s. 

Remark 5. As mentioned previously, in case of the variable J7, there are two alternative ways to 
introduce the boundary condition at x — to the system -formulations (51) and (61), respectively. In 
this way one can construct two alternative dynamic systems of different dimensions (N and A^ — 1). The 
results in Figure 1 and Table 1 show that the stiffness properties of the system corresponding to the 
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Ql/qo = 0.9 


Ql/qo = 0.5 




Vi 


(3) 
<il 




V 


(3; 
1i 


w estimated for the system based on variable w 


(1) 

•km 


6.5c+0 


6.6c+0 


6.8e+0 


1.8e+l 


1.8c+l 


1.8c+l 


( l 2) 
•km 


1.7c+0 


1.7c+0 


1.7c+0 


4.6c+0 


4.7c+0 


4.7c+0 




w estimated for the system based on variable U 


(1) 

•km 


3.0c+0 


3.0e+0 


3.2c+0 


6.0e+0 


6.1e+0 


6.2c+0 


(2) 
•km 


7.5e-l 


7.7e-l 


8.1e-l 


1.5c+0 


1.5e+0 


1.6c+0 




vo estimated for based on condition (51) 


r (i) 

•km 


4.8e+l 


4.8c+l 


4.9e+l 


1.7c+l 


1.7e+l 


1.7e+l 


r m 

•km 


1.2e+l 


1.2c+l 


1.3e+l 


4.3c+0 


4.3e+0 


4.3c+0 




vo estimated for f2( 2 ) based on condition (61) 


T (i) 

•km 


2.3e+l 


2.2c+l 


1.9e+l 


9.6c+0 


1.0e+l 


1.3c+l 


T w 

•km 


5.8c+0 


5.7c+0 


4.7c+0 


2.5e+0 


2.6c+0 


3.5c+0 



Table 1: Values of the parameter from the approximation of the condition ratio (70) for the different 
dynamic systems (67). The computations were provided for e = 1CP 3 

boundary condition (61) are slightly better than those of the system utilizing (51). Indeed, the respective 
parameter w is about two times smaller. One of the possible explanations is the aforementioned difference 
in the systems' sizes: dim A^? = dim a£? + 1 ( see Remark 3). We have checked that the accuracy of 
computations remains practically the same regardless of the system variant. Taking this fact into account, 
we restrict ourselves in the analysis only to the system which employs the boundary condition at point 
x = in the form (61). Thus, from now on all the investigated dynamic systems (for all dependent 
variables) will be based on the same mechanisms for incorporation of the boundary conditions. 

All the results of the stiffness investigation arc collected in the Table 1 and Fig. 1- Fig. 2 The following 
conclusions can be drawn from this data: 

(i) The nonuniform mesh reduces the stiffness approximately up to five times regardless of the solution 
type (Table 1); 

(ii) For the analysed benchmarks it turned out that system stiffness depends on the value of the ratio 
Ql/qo rather than on the particular distribution of the leak-off function (and its behaviour near 
the crack tip); 

(iii) The most important parameter affecting the stiffness properties seems to be the relation between 
the injection flux rate and the leak-off to formation Qi /qg (or equivalently the value of the parameter 
7„ determining the range of the particle velocity variation), as can be clearly seen in Fig. 1; 

(iv) When comparing systems for various dependent variables, the lowest condition ratio gives the 
system built for U (less in order of magnitude in comparison with the others). The worst stiffness 
performance appears for the system corresponding to the il-variable. However this trend becomes 
less distinguishable for the smaller ratios Qi/qo and non-uniform mesh. For example, Fig. lc) 
shows that in some cases f2 may produce lower stiffness than w. 

(v) A value of the regularization parameter e essentially affects the stiffness of dynamic system. This 
is clear when one analyzes Fig. 2, presented only as an example given for a uniform mesh and 
the benchmark qi 1 ' for Qi/qo = 0.9 (see Appendix C). For other combinations of the benchmark 
solutions and different meshes the graphs have similar character. Indeed, the estimation (70) only 
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N 



Figure 2: Effects of different values of e on condition ratio. Type of the benchmark q\ (see Appendix 

C). Condition ratio k = for the w-system for various values of the regularised parameter e in the 

case of uniform grid. There are two behaviours depending on the choice of ./V relative to e. 
a) b) c) 




200 400 600 800 1000 200 400 600 800 1000 200 400 600 800 1000 



TV TV TV 

Figure 3: Maximal relative errors of the solutions computed in different variables 8w, 5U and Sfl, for 
various number of the grid points TV in case of the nonuniform mesh x) n (S = 2). Different values 
of e have been considered. All computations were performed for the benchmark q, 1 ^ for Qi/qo = 0.9. 
Solutions Ui and U n obtained by unitization of the linear and nonlinear regularized conditions (65) and 
(66), respectively 



holds true for sufficiently large TV. The threshold value of TV depends on the chosen e. Thus for 
a fixed number of grid points TV, there is a critical value of the regularised parameter e s (N). By 
taking e < s s (N) one increases the system stiffness appreciably. The relationship between the 
optimal value e s and the number of grid points TV has been roughly estimated from our numerical 
results as s s {N) m g s /N 2 , where the constant g s = pi/' 1 also depended on type of the solution (was 
different for different benchmarks) and on the chosen mesh. 

3.4 Accuracy of the computations 

In this section, we analyse accuracy aspects of numerical computations by MATLAB solver ODE15s 
based on the Rungc-Kutta method specifically dedicated for stiff dynamic systems. This allows us to 
perform a direct comparison between the different approaches (for various dependent variables). The 
predefined accuracy parameters for odel5s, absolute and relative errors of solution, were both set to 
10 -8 which was sufficient to check the maximal achievable accuracy for different dynamic systems. Note 
that the stiffness properties of the systems are not directly related to the accuracy of computations: the 
best stiffness performance of a solver does not immediately correspond to the best accuracy of solution. 
As will be shown later, both concurrent tendencies can be observed. 
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Comparison of conditions (58), (65), (66) 
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= 0.9 




Ol Inn 


= 0.5 






IQ-' Z 


10~ 4 


icr b 


10~ 2 


io- 4 


10~ b 


6U* 




1.6c-l 


1.4e-l 


1.3e-l 


6.1c-3 


3.7e-3 


3.7c-3 


T w 


1.4c-l 


7.6e-2 


6.3e-2 


4.5c-3 


9.3e-5 


8.9c-5 


SUi 


r (i) 


5.0e-2 


1.4e-2 


1.7c- 2 


1.2c-5 


l.lc-5 


l.lc-5 


T w 


5.0e-2 


1.7c-3 


2.0e-3 


4.9e-5 


8.2e-5 


8.7c-5 


8U n 


T U) 


4.4e-2 


1.2e-2 


1.3c-2 


2.2e-5 


1.3e-5 


1.3c-5 


r w 


4.4e-2 


9.9e-4 


1.8e-3 


4.2e-5 


8.2e-5 


8.7e-5 



Table 2: Comparison of the accuracy of the solution of dynamical system based on variable U. The 
results depicted by [/* refer to the regularized boundary condition based on one asymptotic term, while 
those denoted by Ui and U n correspond to two terms approximation (linear (65) and nonlinear (66), 

(2) 

respectively). Other problem parameters: N = 100, 5 = 2 for the mesh x m . 

Before we compare different approaches in terms of their accuracy, let us recall two alternative ways 
to define the regularized boundary condition at the end point x = 1 — s. The first one is based on the 
e-regularization technique, as it was defined in Linkov (2011d) (sec also Mishuris ct al. (2012)). The 
condition is equivalent to implementing the test function into the crack tip element taken in the form 
of a single (leading) term expansion of the solution asymptotics near the crack tip. For the variable 
U this condition is defined by (58). The second approach to formulate the regularized condition is to 
take into account the first two terms of asymptotics as described in section 3.1 (compare equations: (62) 
and (63)) 8 . Finally, in the case of U, this condition may be implemented in the non-linear form (66). 
One can expect that the two terms conditions would have a clear advantage, at least in cases when the 
solution smoothness near the crack tip is deteriorated due to the singularity of the leak-off function. 

The results of the computations presented in Table 2 confirm such prediction. We compare only 
conditions for U, as originally the e-regularization technique was introduced for this variable. Indeed, 
the relative errors of the solutions 5Ui or SU n are at least one order of magnitude lower than that in 
the case of 5U*, corresponding to the formulation based on (58). 9 Surprisingly, for the variants of the 
non-singular leak-off function, the improvement is even more pronounced (especially for uniform mesh). 
Thus, the implementation of further asymptotic terms in the formulation of the regularized boundary 
condition always leads to better results, as one could predict. 

We also made the computations for three different benchmarks reported in Mishuris et al. (2012). 
They correspond to the leak-off function vanishing near the crack tip. It turned out that computa- 
tional error corresponding to the modified form of the regularized conditions (based on two terms of 
asymptotics) was always two orders of magnitude lower than that reported in the previous paper. 

On the other hand, there is no difference observed between the solutions SUi but 5U n at least for 
those two benchmarks and the choice of the parameters (N = 100). However, as we will show later, for 
large numbers of nodal points, or more severe leak-off function relationship (Qi/qo ~ 1), the nonlinear 
formulation of the condition clearly manifests its advantage. 

From now on only the regularized conditions based on two asymptotic terms (63) will be utilized. 
Additionally, for variable U , two different forms, linear (65) and nonlinear (66), will be adopted. For w 
and Q two formulations, (59) and (60), which are equivalent to the condition (58) could not compete 
with their more accurate analogue (63) in terms of solution accuracy and will not be considered. 

8 One can argue that such approach is equivalent to implementation of the multiscale asymptotics in spirit of Garagash 
et al. (2011) 

9 Here and everywhere later, by Sf we understand the maximal value of the relative error of the function / over all 
discretized independent variables (<5/ = ||<5/||oo). 
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Dynamic system built on the variable w 
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Ql/qo = 0.5 



Table 3: Performance of the solver based on the dependent variable w for number of nodal points 
N = 100 and various benchmarks. Values of the regularized parameter are e = 5 • 10~ 3 and e = 10~ 3 
for the meshes for x£} and x$ , respectively 



Graphs presented in Fig. 3 illustrate some peculiarities of the computational process. Here, the 
maximal relative errors of the solutions (over the time and space) as functions of the number of mesh 
points N (Sf = 6f(N)), are presented for different variables. The considered benchmark assumes 
Ql/qo = 0.9 (see Appendix C). Three different values of the regularization parameter e = 10 -3 , 10~ 4 
and 10~ 5 were chosen. 

In Fig. 3 two basic tendencies can be observed. The first one is the monotonous error decrease with 

growing N, up to some stabilization level. This level is different for different dependent variables and 

values of e, and in some cases is reached for N > 1000 (and thus cannot be identified in the figure). The 

second trend is discernible when comparing results for different values of e. Namely, it turns out that 

for each dependent variable there exists an optimal e minimizing the solution error. This value however 

depends on N. One can approximate the optimal value of the parameter in a way similar to that for the 

condition ratio: e a (N) rj g a /N 2 . The constant g a = gt/^ does not coincide (in fact may differ in orders) 
( f ) 

with g s found earlier in the stiffness analysis. It is not a surprise that the optimal stiffness properties 
and the maximal solution accuracy are not achieved for the same values of the regularization parameter. 
Indeed to increase computational accuracy one needs to decrease e and increase number of the nodal 
points N. However, both of these leads to increase of the condition ratio. 

Note that the relative errors of respective dependent variables cannot be compared directly. Indeed, 
even if the errors for w and U are interrelated via the evident relationship SU = 3dw, their comparison 
with SQ necessitates an additional postprocessing of the latter. This process, in turn, may introduce its 
own error. On the other hand, there exists a common component of the solutions, the crack length SL, 
which can be naturally used for such comparison. 

Below we adopt the following strategy for performance test for different dynamic systems. First, 
we set the number of nodal points, N, to 100. Next, for each of the dependent variables we accept 
optimal (for N — 100) values of the regularization parameter e. It turned out that the optimal e differs 
slightly depending on the type of mesh chosen and the benchmark variant. The general trend for e can 
be identified for different meshes (for Xm it is always smaller than for ). However, the sensitivity 
to the benchmark type is low. The results of computations described by various accuracy measures are 
collected in Tables 3-6 (the optimal values of e are specified in the captions). We present there: the 
relative error of solution Sf, the absolute error of solution Af and the relative error of the crack length 
SL. The following conclusions can be drawn from this data: 

(i) There is a trend of simultaneous increase of the ratio Qi/qo and the relative errors of dependent 
variables Sf. However this tendency is not in place (or may be even reversed) when analyzing SL. 
The solution accuracy is affected much stronger by the value of Qi/qo itself, than by the leak-off 
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Table 4: Numerical results for the system built on the dependent variable U with the linear regularized 
condition (65) for N = 100 and different e for the uniform and nonuniform meshes (e = 10~ 4 and 
e = 10 -5 , respectively). 
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Table 5: Results for the solver based on the dependent variable U for nonlinear regularized condition 
for N = 100 and various benchmarks. Values of the regularized parameter are e = 10~ 4 and e = 10~ 5 
for the meshes for i„ and , respectively 
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Table 6: Computation accuracy for the solver based on the dependent variable ft for N = 100 and 
e = 10~ 2 for and e = 5 • 10~ 3 for x$. 
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function behaviour near the crack tip. 

(ii) In case of dependent variable U, the way in which the regularized boundary condition is intro- 
duced (linear or non-linear) does not play an essential role for the benchmarks and ranges of the 
parameters under consideration in the Tables 3-6. However, there are exceptions to this rule. 
One of them can be seen in Fig. 3c) for large values of N, where the non-linear condition proves 
its superiority. Another case will be presented in the end of this section. 

(iii) When comparing the common accuracy parameter SL, the dynamic system for fi gives the best 
results. Dynamic system for w is the worst performing scheme and comparable to the one for U 
in a few cases. 

(iv) Since f2 vanishes near the crack tip faster than others variables, one could expect the worst relative 
error in this case. Surprisingly, even when contrasting the relative (incomparable) errors of the 
respective dependent variables with each other, the system for f2 seems to be the best choice. The 
advantage of D, over w and U is especially pronounced for the benchmarks variants with a higher 
ratio Qi/qo- The worst option from the accuracy point of view is the system based on w. 

(v) Better solution accuracy is obtained for the non-uniform mesh in almost every case. 

We have not observed any significant difference between the time step strategies chosen by the 
ODE15s solver for different dynamic systems. The number of steps and the main trends were similar. 
For this reason we have not presented any details in the tables. 

To visualize the results reported in the Tables 3-6, and to complement those presented in Fig. 3, in 
Fig. 4 we show the relative errors of the crack length SL computed by different dynamic systems (built on 
different variables). The same benchmark and the values of all other problem parameters as previously 
discussed in Fig. 3 were considered. If the trends for the different relative errors of the solution Sf and 
the crack length SL in case of e = 10 -3 look similar, the results for the smaller value of the parameter 
are rather surprising. Indeed, the relative errors SUi and SU n are smaller than Sw and (557, while the 
error SL computed for U no longer follows this trend. 

This paradox needs an explanation. A trivial one could be that the maximal error for the solution 
(Sf) is not situated near the crack tip but inside the computational domain. To verify this hypothesis 
and to give a prospective reader a clear picture of the distribution of the solution error in time and space, 
we present, in Fig. 5 and Fig. 6, the corresponding absolute and relative errors computed for nonuniform 
mesh built on N = 100 nodal points with the corresponding optimal regularized parameters discussed 
after Fig. 3. As it follows from Fig. 6, the maximum of the relative error is always achieved near the 
crack tip (maxt 8f(t, 1 — e))). Hence, the initial guess has not been confirmed. On the other hand, the 
values of the relative errors Sf and the respective SL^ are directly interrelated. The following analysis 
identifies this relationship. 



Using (62), after some algebra, one has the estimate: 



(/) 

Sf a Se[ f) + (Self - Se[ f) )%, e a2 ~ ai . (71) 



.(/) 

For the benchmark and Qi/qo = 0.9 (see Appendix C) which always provides the worst accuracy in 
our computations, one can conclude 

Swk isu^Sw + ^(5w 1 -Svo )^e, (72) 
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Figure 4: Distribution of relative errors of the fracture length SL computed by solvers based on different 
depended variables (w, U and il). When dependent variable U is considered, two different regularised 
boundary conditions arc in use: (65) for Ui and (66) for U n . Other parameters are the same as in Fig 3. 
Zoom picture within the Fig. 4 b) corresponds to the sharp minimum of SL for the variable Ui. 
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Figure 5: Absolute error for solutions w, U n and Q computed for benchmark q^ with ratio Qi/qi = 0.9 
and nonuniform mesh (S = 2) with N = 100 nodal points. Other parameters: e = 10~ 3 for w, 
e = 5 • 10~ 3 for £1, and e = 10~ 5 for U n . 



5n*>6w + — (Swi-Sw )^. (73) 

Finally, from (31) we can derive 

SL « ^Sw . (74) 

The last relationship has also been verified numerically by evaluating the values of the constant wq in the 
postprocessing procedure using the computed solution (w, U or Q) and the corresponding regularized 
boundary condition (compare (62) and (63)). 

It is clear from relations (72) - (73) that the relative errors of the respective dependent variables also 
depend on the quality of approximation of the second term in the regularized boundary condition (62). 
This explains the surprising relationship between SL and the respective 5f. 

Interestingly, the results presented in Fig. 4 show that the value of e which provides the lowest relative 
error, Sf, of the dependent variable / does not necessarily give the best accuracy of the crack length 
SL. Moreover, the relation e = £l(N) is much more sensitive to the variation of N than e = Ef(N). 
Indeed, one can observe sharp minima (see Fig. 4 a) and b)) while there is no such phenomenon in the 
respective graphs for Sf (see Fig. 3). To demonstrate that the peaks are not computational artifacts, 
we also include a small zoom of the corresponding area of the figure Fig. 4 b). 

To complete the accuracy analysis, let us consider some critical regime of crack propagation. Namely, 
assume that the leak-off flux almost entirely balances the volume of fluid injected into the crack. Indeed, 
when taking the Carter type benchmark (92) b\ = 62 = 1, one obtains the fluid balance ratio Qi/qo = 
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Figure 6: Relative error of the solutions w, U n and Q computed on the corresponding solvers for the 
same parameters as in Fig. 5. 
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Table 7: Accuracy parameters for the limiting (critical) variant of the benchmark solution (Qi/qo = 
0.9857, 7„ = 2.07) computed for different meshes composed of TV = 100 nodal points. The blank positions 
in the table correspond to the case when the solver ODE15s could not complete the computations in a 
reasonable time. 

0.9857. This gives a very strong variation of the particle velocity function along the crack length 
(7, = 2.07). 

In view of the previous conclusion on the ratio Qi / qo influence on solution accuracy (which in fact 
confirms the observations from Mishuris ct al. (2012)), one can predict that the solution error will increase 
appreciably in comparison with the figures shown in Tables 3-6. In order to verify this assertion the 
computations were made for respective dynamic systems (the system for U was analyzed again for two 
forms of the regularized boundary condition). Both types of meshes, the uniform and the non- uniform, 
were utilized, each composed of 100 nodal points (N = 100). Four different values of the regularization 
parameter e , ranging from 10~ 5 to 10 -2 , were analyzed. The results of the computations described by 
respective accuracy parameters are presented in Table 7. Here, the symbols SUi and SU n stand for the 
relative error of U obtained for the conditions (65) and (66), respectively. The subscript of 5L informs 
us which dynamic system the corresponding result was obtained for. 

The data in the table shows that, regardless of the accuracy parameter, the solution error increased 
at least one order of magnitude, as compared to the values from Tables 3-6. The lowest deterioration 
of the solution accuracy was obtained for fi, which proves the best overall performance of the system 
built for this variable. Especially impressive is its advantage when comparing the errors of crack length 
estimation. In all considered cases SLq is at least two orders of magnitude lower than SL for other 
dependent variables. 

In this critical variant of benchmark solution, the non-linear regularized boundary condition (66) for 
U gives, in most cases, much better performance than its linear counterpart (65) (compare with the 
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discussion after the Tables 2 and 3-6). Finally, the non- uniform mesh seems to be a better choice from 
the point of view of accuracy. 

4 Conclusions 

In this paper, we have revisited the PKN model of hydrofracturing providing a comprehensive overview 
of the known results and supplementing the state of the art with: (i) analysis of the leak-off term in the 
Reynolds equation in its normalised form for the case of Carter model; (ii) complete asymptotic analysis 
of the solution for this model; (hi) new dependent variable, f2, and the resulting problem reformulation; 
(iv) new form of the regularized boundary condition. 

The main outcome of the paper is a thorough analysis of the dynamic systems constructed on different 
dependent and independent variables in terms of the system stiffness, accuracy and efficiency of the 
computations. 

The results show that the approach proposed in Linkov (2011d); Mishuris ct al. (2012) may be efficient 
even in the cases when the smoothness of the particle velocity near the crack tip is disturbed by the 
singular leak-off function. However, to guarantee the best performance of the solver, a modified way to 
impose the regularized boundary condition is advisable. Moreover, this modified formulation proves to 
be beneficial also for the non-singular leak-off regimes. 

The stiffness analysis performed shows that the dynamic system stiffness can be reduced by a proper 
manner of spatial meshing (proper choice of independent variable). Simultaneously, the non- uniform 
mesh may appreciably reduce computational error. As could be expected, the mesh optimal for the 
system stiffness is not the best choice from the accuracy point of view. Thus, a compromise between the 
efficiency and accuracy of computations should be made, when accepting the choice of spatial meshing. 

The new independent variable, f2, shows its superiority over other formulations. Indeed, it consid- 
erably improves the accuracy of computations (the fracture volume, Q(t,x) and especially the fracture 
length, L(t)). The dynamic system built for this variable exhibits much less sensitivity to the process 
parameters than the respective systems for w and U. However, its stiffness may be slightly increased in 
comparison with the alternative variants (for U and w) . 

Although, the presented analysis concerns the 1-D PKN formulation, at least some of the findings 
may be utilised in the more advanced cases. 
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Appendices 

A Carter's leak-off function in the normalised formulation 

Consider the transformation of the Carter law described by (4) when applying the normalization (18). 
Assume that: 

J— = ^k+R(t,x), (75) 

yft — T(X) V 1 — X 

where function D(t) is defined in (22) while the remainder R is estimated later in (78). 
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To find function D(t), and thus to obtain an exact form of equation (22), it is enough to compute 
the limit 

D 2 {t) = lim fr. (76) 

x-¥l t — T\X) 

This can be done by utilising L'Hopital's rule with taking into account that x — > L{t) as x — > 1, 

t(x) = t (L(t)x) = L _1 (£(i):r), (77) 

and that the crack length is a smooth function of time (L <G C 1 at least). The last fact immediately 
follows from the problem formulation in terms of evolution system (32). 

Having the value of D(t) we can estimate the remainder R(t, x) when x — > 1, or, what it is equivalent 
to when x —> l(t) (or t — > t(x)). For this reason, we search for a parameter £ ^ which guarantees that 
the limit 

A=Hm *M) =lim i ( D(t) L(ty ix) 



s->i(l-s)« 2£(1 -i)?- 1 V(l-i) 3/2 (t-r(a;)) 3 / 2 

does not turn to zero or infinity. Due to this assumption, we can write 

-=J== = -^L + A(l -xf + o ((1 - , (78) 

^/t - t(x) Vl — 35 

when x — > 1, or equivalently a; — > Z(t). Taking the last estimate into account A can be expressed as: 

. 1 / D{t) L(ty(x) D(t) \ AL(t) (l-5)r / (x) /i ,^ 
^ = lim t-tt, tt^-t t -ttt^ , / — f==== -r^ lim ^ '-. \ ' (l + oil )). 

Now, on substitution of r'(a;) = 1/ L'(t) at x = L(t) and (76) into the limit one has: 

A = V ° {t) ( 1 WM ^l ^(t)J 2 W 

s->i 2^(1 - 2)C-V2 ^1-5 t-r(x) J 2£L'(t) ' 

Applying (76) and (22) here gives: 



1 + 2^ = £>(f) / _J^_ _ £(t)r'(z) J>(t) \ _ A£>(W) T ; (a;)VT~~x 

2£ s ™ 2£(1 - s)€-Va ^ ! _ j ^ _ T ( a; ) ^T^i J 2£ s™ _ T (x) ' 

By repeating the same process one more time we have: 

(2 + 20 A = lim *W f 1 - WW)) . 

Finally by eliminating the square root with use of (78) we obtain (after some algebra) 

I- L(t)T'(x)D 2 (t) 

This relationship gives a finite value of ^4 if and only if £ = 1/2 and, as a result, we find: 



A = -D 3 (t)L 2 (t) T "(L(t)) = --^M J^Q. 
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B Asymptotics of the solutions for different leak-off functions 

Asymptotic expansion for the crack opening and the fluid velocity near the crack tip in the normalised 
variables (18) can be written in the following general forms: 

N 

w(t,x) = y ^w j {t){l - x) a * + <3((1 -x) e ™), (79) 

3=0 

and 

N 

V(t,x) = J2v j (t)(l-x)h +O((l-a0"), x^l, (80) 

3=0 

with g w > a n , qv > /3 n , ao = 1/3, fio = and some increasing sequences oto,a.\, . . . ,a n and 
Po, Pi, . . . , (3 n . Note that the asymptotics are related to each other by the speed equation (19) and 
thus, regardless of the chosen leak-off function, we can write 

N N N N 

Y^Vim-xfi + ... = E E^^CK^Hwa-a)"^^- 1 . (»i) 

In line with the discussion after equation (16), we are interested only in the terms such that f3j < 1, 
restricting ourselves to the smallest gy > 1, since the values of j3j are combinations of a sum of three 
consequent components of the exponents ay. However, since «o is known (cto = 1/3), one can write 
(compare with (17)): 

v ^ - m w3M (82) 

Vi(t) = j^(ai + l)w*(t)w 1 (t), /8i=ai-i. (83) 

To continue the process one now needs to compute the value of the exponent a\ as it is not clear a 
priori which value determining the next exponent /3 2 = min{2/3 + a 2 , 1/3 + 2a!} is larger. To do so let 
us rewrite the continuity equation (20) in the form: 

dw V a (t) 8w 1 d{w(V - V)) 

X) ^ = W) dx qi(t ' x) - (84) 

Here, the terms on the left-hand side of the equation are always bounded near the crack tip, while those 
on the right-hand side can behave differently depending on the chosen leak-off function. 
Consider the following three cases of qi behaviour. 
(i) Assume first that 

qi(t,x) = o(w(t,x)) , x 1. 

This case naturally includes the impermeable rock formation. Analysing the leading order terms in the 
equation (84), it is clear that w(V - V) = 0((1 - x) 4/3 ), as x -> 1. This, in turn, is only possible for 
pi = 1 and, therefore, «i = 4/3. Finally, comparing the left-hand side and the right-hand side of the 
equation we obtain: 

<(t) = ^M(Vo(t)+4V 1 (t)), V 1 (t) = ^w 2 (t)w 1 (t). (85) 
This case has been considered in Linkov (2011d) and Mishuris et al. (2012). 
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(ii) If we assume that the leak-off function is estimated by the solution as 0(w(t, x)) , or equivalcntly; 

qi(t,x)~r(t)wo(t)(l-x)V s , x->l, 

then the previous results related to the values of u\ and fa and, therefore, the equation (85)2 remain 
the same, while the first one changes to 

<(t) - 3^«>o(t)fa(t) + 4Vi(t)) - T(t)w (t). (86) 

This case corresponds to (21)3 when C32 = and T(t) = kC^\{t). 
(m) The leak-off function in a general form: 

qi{t,x) = $(t)(l-z) e +o((l-z) 1/3 ), x -> 1, 

where — 1/2 < 6> < 1/3. Here, one can conclude that iy(Vo — V) = 0((1 — x) 1+e ). as a; — s- 1 or cquivalcntly, 
fa = 6 + 2/3, and ai = 1 + 0. Moreover, in this case: 

(1 + 6)w Q V 1 = L(t)$(t), V^t) = J—(e+f) wl(t)w!{t), (87) 



L(t) V 3 

and, thus 

_ 3L 2 {t)$(t) 
Wl{t) ~ {4 + 39)(l + 6)w 3 Q (t)- 

Note, that as one would expect, the particle velocity function is not smooth in this case near the 
crack tip, its derivative is unbounded and exhibits the following behaviour: 

d ^ = o{{i- x f-y% 

To formulate the equation similar to (85) 1 or (86), one needs to continue asymptotic analysis of the 
equation (84) incorporating the available information. Apart from the fact that the analysis can be done 
in the general case, we restrict ourselves only to three variants used from the beginning (compare (4)), 
respectively: 9 = 0, 6 = 1/3 - 1/2 = -1/6 and 6 = -1/2. 

When 6 = 0, ol\ = 1 and fli = 2/3, returning to the equation (81), one concludes that /?2 > 1 and, 
therefore, 

1 

MAt) 



W 'S) = 7rr7^Mt)v G {t). (89) 



This case corresponds to (21)3 when $(£) = C3 (t)wo(t) and C3 = 0. 

If 9 = —1/6, then a.\ = 5/6 and /3i = 1/2. In this case the function $(i) can be written as 
$(t) = C2D(t)wo(t) (compare to (21)2) and again equation (81) gives P2 > T while equation (84) leads 
to 

<(t) = 3Mj( w °(tWo(t) +4w 1 (t)V 1 (t)). (90) 

Summarizing, in both mentioned above cases, there exists a single term in asymptotics of the particle 
velocity which has singular derivative near the crack tip. Moreover, those terms (w\ and V\, respectively) 
are fully defined by the leak-off function $(t) and the coefficient wq in front of the leading term for the 
crack opening in (88) and (87) 1. 

The situation changes dramatically when 6 = —1/2 (Carter law). We now have a.\ = 1/2 and 
pi = 1/6 and <I>(i) = C\D(t). In this case, however, fa < 1 and we need to continue the asymptotic 
analysis further to evaluate all terms of the particle velocity which exhibit non-smooth behaviour near 
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the crack tip. We omit the details of the derivation, presenting only the final result in a compact form. 
The first six exponents in the asymptotic expansions (79) and (80), that introduce the singularity of w x , 
are: 

a j = o + «' = i = 1 > 2 .---»6, 



2 6 J 6 



and 



„ _ 12 _ _270 _ 9768 2097252 _ 1081254096 

7 ) *>2 49 j ^3 ~ 343 j 120 05 ' ^ 5 — 924385 ' 

24 „/, _ 828 „/, _ 5136 „/, _ 1234512 



^1=2, lfe = -f, V^3 = ^f, ^4 = "^, ^5- 

C Benchmark Solutions 

There are several benchmarks in the literature to be utilized for investigation of the numerical algorithms. 
Benchmark solutions for impermeable rock have been constructed in Kemp (1989); Linkov (2011d), while 
that corresponding to the non-zero leak-off model with qi vanishing at a crack tip has been analyzed in 
Mishuris et al. (2012). 

In this paper, we introduce three different analytical benchmark solutions corresponding to the 
representations (21). Moreover, for each of the leak-off functions under consideration we take two 
different relationships between the injection flux rate qo and the leak-off to formation qi . In this way six 
different benchmark solutions are analyzed. 

In order to formulate the benchmark solutions let us assume the following form of the crack opening 
function: 

w(t,x) = W (l + tyh(x), W = y|(3 7 +l), (91) 

where 7 is an arbitrary parameter, and the function h(x) (0 < x < 1) is given by: 

h(x) = (1-x)* +b 1 (l-x) Xl +b 2 (l~x) X2 . (92) 

The choice of the next powers 1/3 < Ai < A 2 will depend on the leak-off variant from (4). On consecutive 
substitutions of (91)-(92) into the relations (19), (24), (29) and (31) one obtains the remaining benchmark 
quantities: 

L(t) = (l + t)^, V (t,x) = -W^l + t) 3 - 2 T 1 h 2 (x)^. (93) 



qo{t ) = -W i (l+t)^(h 3 ^) \ x=0 . 



(94) 



It can be easily checked that for Ai = 1/2 and A 2 = 4/3 the leak-off function incorporates a square 
root singular term of type (21)i. By setting Ai = 5/6 and A 2 = 4/3 we comply with representation (21) 2 . 
Although in both of these cases <7*( 2 ) exhibits a singular behaviour at the crack tip, it does not detract 
from the applicability of our benchmarks. Finally, when using Ai = 4/3 and A 2 = 7/3, the benchmark 
gives a non-singular leak-off function in the form (21)3. 

Note also, that by manipulating with the value of 7 one can simulate some very specific regimes of 
crack propagation. For example 7 = 1/5 corresponds to the constant injection flux rate, while 7=1/3 
gives a constant crack propagation speed. For our computations we always set the value of 7 = 1/5. 
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a) b) c) 




x x 



Figure 7: Distributions of the leak-off functions qi(t,x) and the respective particle velocity V(t,x) over 
x G (0, 1) at initial time t = 0. 



Choosing appropriate values b\ and bi one can change the relation between the amount of fluid loss 
to formation and the injection rate. This ratio can be defined by the measure, Qi/qo, where Qi is the 
total volume of leak-off J Q qidx. It is important to note that this measure decreases in time, from its 
maximum value to zero, for all chosen benchmarks. Thus, taking the maximal value high enough and 
tracing the solution accuracy in time, one can analyse performance of the algorithm for any possible 
value of the parameter. We consider two variants of Qi/qo, one where fluid injection doubles the size 
of total fluid loss, and a second where the total fluid loss is close to injection rate. The values of the 
corresponding parameters b\, &2 are presented in Table 8. 





Qi/qo = 0-9 


Qi/qo = 0.5 




1l 


(3) 








bx 


0.1 


0.19 


0.74 


0.02 


0.03 


0.15 


b 2 


0.5 


0.41 


-0.13 


0.1 


0.08 


-0.02 


Jv 


1.69 


1.70 


1.65 


0.55 


0.56 


0.55 



Table 8: The values of parameters b\ and 62 for different benchmark solutions modelling desired leak-off 
to fluid injection ratios. 



Additionally one can compute a parameter 7„ defined in Mishuris et al. (2012) as a measure of the 
uniformity of fluid velocity distribution: 



7„ = [max(V(i, x)) — mm(V(t, x))] 



, -1 







Interestingly, this measure is directly correlated with the leak-off ratio Qi/qo- 

In Fig. 7 the distributions of the leak-off functions and the corresponding particle velocities for the 
respective benchmarks are presented. It shows that the velocity near the crack tip depends strongly on 
the benchmark variant. To highlight this fact, a zoom picture is placed in the Fig. 7 b). 

Note that the benchmark qj 1 ' is worse, in a sense, than the original Carter's model as it contains 
additional singular terms of the leak-off function. These terms are absent in the normalised Carter's law 
as it follows from Appendix B. 
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